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ABSTRACT 

We present the implementation of a three-dimensional, second order accurate 
Godunov-type algorithm for magneto- hydrodynamic (MHD), in the adaptive- 
mesh-refinement (AMR) cosmological code CHARM. The algorithm is based on 
the full 12-solve spatially unsplit Corner- Transport- Upwind (CTU) scheme. The 
fluid quantities are cell-centered and are updated using the Piecewise-Parabolic- 
Method (PPM), while the magnetic field variables are face-centered and are 
evolved through application of the Stokes theorem on cell edges via a Constrained- 
Transport (CT) method. The multidimensional MHD source terms required in 
the predictor step for high-order accuracy are applied in a simplified form which 
reduces their complexity in three dimensions without loss of accuracy or robust- 
ness. The algorithm is implemented on an AMR framework which requires spe- 
cific synchronization steps across refinement levels. These include face-centered 
restriction and prolongation operations and a reflux- curl operation, which main- 
tains a solenoidal magnetic field across refinement boundaries. The code is tested 
against a large suite of test problems, including convergence tests in smooth 
flows, shock-tube tests, classical two- and three-dimensional MHD tests, a three- 
dimensional shock-cloud interaction problem and the formation of a cluster of 
galaxies in a fully cosmological context. The magnetic field divergence is shown 
to remain negligible throughout. 

Subject headings: cosmology: theory — methods: numerical — magnetohydro- 
dynamics: MHD 
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Introduction 



Magnetic fields are a common feature of cosmic plasmas, from the interplanetary medium 
and the atmospheres of stars, to the interstellar medium of galaxies and the baryonic 
gas in the largest structures of the universe such as clust e rs and voids of galaxi es (e.g., 
Zel'dovich erZlllQSsi ISernet et allboosi Iciarke et allboOll : iNeronov fc Vovklboiol ). 



Their origin is discussed in several papers and diffe rent processes are lik ely responsible 
for the magnetization of different environments (see, e.g. iMiniati fc Bellll201ll and references 
therein). In general weak rotational electrostatic fields are required, which are normally sup- 
pressed by the high conductivity of astrophysical plasmas, but which can nevertheless arise 
under special conditions. The magnetic field can then evolve considerably and amplifica- 
tion of an initially weak seed by many orders of ma gnitude is plausible, particularly when 



the flow is highly turbulent (jZel'dovich et al.lll983[ ). Magnetic fields affect the dynamics 



of a system directly through the Lorentz force and indirectly through their impact on the 
plasm a microscopic properties, e.g. the thermal conductivity or el ectric resistivity ( ISpitzer 
19651 ). or the transport of high energy particles ( ISchlickeiserl 120021 ) . both with conspicuous 
macroscopic effects. Thus magnetic fields are a crucial component of astrophysical plasmas 
although perhaps due to the complexity they introduce, progress in characterizing their role 
has been relatively slow, particularly in certain areas of astrophysics. Due to such complex- 
ities, particularly during highly nonlinear regimes, accurate and efficient numerical methods 
are valuable for studying the evolution of magnetized systems. 

A simplified description of a magnetized fluid is provided by the equations of ideal 
magneto-hydrodynamics (MHD). This approximation is valid when the following hierarchy 
of scales is satisfied: f^/p ^ A ^ L, where l^fp is the mean-free-path of the fluid p articles, A 
is the characteristic scale of the problem of interest and L is the size of the system (IGinzburg 
19791 ). The mfp can be suppressed considerably in the direction perpendicular to the mag- 
netic field lines, but in parallel directions one relies typically on collisions to guarantee the 
fluid approximation. In reality microscopic plasma processes or a tangled component of the 
magnetic field can also provide the required viscosity. Notably the MHD approximation 
prevents the occurrence of kinetic processes which, however, can sometimes by represented 
by source terms on the RHS of the MHD equations. 



When dissipative terms can be neglected, the MHD equations are in ideal form and 
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read (iLandau &: Lifshitall984j ) 



dt 
dpe 



dp dpuj 
dt dxj 



d 



+ (pUiUj + pSij - BjBi) = 0, 
d 



dt dxj 



0. 



(2) 
(3) 



Here p is the gas density, Ui and Bi the components of the velocity and magnetic field 
vectors, respectively, p = Pg + \B'^ is the total sum of the gas and magnetic pressures, 
e = + Gth + \— is the total specific energy density, 5ij Kronecker's delta and summation 
over repeated indices is assumed. The thermal energy is related to the pressure through 
a 7-law equation of state, eth = Pgl pipi — !)• The magnetic field evolution is described by 
Faraday's equation, with the electric field given by Ohm's law. In the limit of negligible 

ds are those induced by motions of the magnetized fiuid and 



resistivity the only electric fie^ 
the induction equation reads (ILandau fc Lifshitall984f ) 



dB, 



dEk 



d 



—e 



dt dxj dxj 

where Eijk is the fully antisymmetric tensor of rank 3 and £012 



{ujBi 



BjUi 



(4) 



1. 



The resistive terms neglected in Eq. which are responsi ble for the diffusion 

of the magnetic field, can be readily recovered (ISamtaney et al.ll2005[ ). although for most 
purposes their neglect is safe in astrophysical plasmas. 

There are various approaches to solve numerically the equations of MHD. The one we 
follow in this paper is based on the extension to MHD of conservative methods for hyperbolic 
systems of equations, particularly Godunov's methods for hydrodynamics. This approach is 
met with two difficulties, however. First, care much be taken because the system of MHD 
equations is not strictly hyperbolic. Th is can be dealt w ith by "renormalizing" the eigen- 
values and eigenvectors of the system (IBrio &: Wul 119881 ). Second, and most importantly 
the solenoidal constraint, V • B = 0, must be enforced or else the solution will contain 



artifacts (IBrackbill fc BarnesI Il980l ) . Unfortunately numerical schemes designed for pure 
hydrodynamics do not fulfill such requirement and the above constraint must be enforced 
separately. Different ways to do so have been proposed. In one approach the magnetic 
field is defined together with all other fiuid variables at cell centers and the non-solenoidal 
componen t is removed through a Hodge- Helmholtz pro j ection method, typically once per 
time-step (IBrackbill &: BarnesI Il980t IZachary et al.lll994j : iRyu et al.lll995l ). In a variation of 
this approach the projection operation is performed on the magne tic field variables ext rap- 
olated to the cell faces which are used to define the MHD fluxes (ICrockett et al.ll2005[ ). It 
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is argued in ICrockett et al.l ( 120051 ) that this type of projection is mathematically more con- 
sistent with the solenoidal requirement because it is the fluxes defined on the cell faces that 
get differentiated to compute the flux updates. In a second approach, the MHD equations 
are cast in a special 8-wave non- conservative formulation. This allows for the non-solenoidal 
component of the magnetic fleld to be explicitly trac ked and suitably damped as it is advected 
by the flow JPowell et al.lll999l : IPedner et al1l2002h . 



Finally in th e Constrain ed Transport (CT) approach, the discretization strategy origi- 



nally proposed by lYed (119661 ) in the context of Maxwell equations is used, in which the mag- 



netic fleld is deflned at face centers, the electric fleld used to update the latter is deflned at 
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this approach the rate of change of the magnetic flux at cell faces is given by the circulation 
of the electric fleld along the cell edges which deflne the boundary of the corresponding 
face. Thus the solenoidal character o f the magnetic hel d is ensured by Stokes' theorem 
down to machine precision. Recently [Gardiner fc Stond (120051 120081 ) have developed an 
unsplit version of the CT algorithm, extending to the case of MHD the Co rner Tra r isport 
Up wind (^CTU ) meth od for directionally unsplit hydrodynamics proposed in IColellal (Il990l ) 
and ISaltzmanI (119941 ). The use of a directionally unsplit algorithm has proven quite at- 
tractive, particularly because it appears to be better suited fo r modeling turbulent flows 
(Bell et al) and in the presence of source terms ( lLevequall998l ). Most importantly, due to 
the solenoidal constraint, the MHD equations contain intrinsically multidimensional terms 
which require a directionall y unsplit formulation if o ne is to achieve high order accuracy in 
multidimensional problems ( iGardiner fc Stond l2005l ). Thus the extension of CTU to MHD 



has also attracted considerable interest in the astrophysical community (^e. g. iTeyssier et al. 
20061 : iFromang et"aDl2006l : llee fc Dean^bood : iMignone fc TzeferacoslEoiol ). 



In this paper we describ e the implementati o n of a versi on of the CTU + C T scheme that 
closely resembles the one of iGardiner &: Stond (120081 ) and IStone et al.l (|2008[ ) . Our scheme 
differs from theirs, however, in two respects. First we have chosen to use the full 12-solve CTU 
scheme instead of the simpler 6-solve scheme. The reason for this choice is due to the larger 
CFL number that the full CTU sch eme can afford (CFL=1), as compared to the simplifled 
version (CFL=0.5). As indicated by IGardiner fc Stond (120081 ) the computational cost of the 
two versions of the CTU scheme is roughly the same for pure MHD calculations, because 
the factor of two fewer Riemann solvers comes with twice as many steps to achieve the same 
solution time. However, for multi-physics applications that we have in mind other expensive 
solve rs are executed each time step, whose cost grows with the number of time-steps. Note 
that iTeyssier et al.l (|2006[ ) use also a simpler version of the CTU scheme, although in their 
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case a slightly less restrictive condition on the time-step is nominally allowed, i.e. CFL< 0.7. 
On the other hand, Lee & Deane (2009) have introduced a different approach for computing 
the transverse flux updates of the fluid variables in their directionally unsplit method, which 
relies on characteristic tracing alone and does not require intermediate Riemann solvers 
(except for the magnetic field intermediate updates). For this reason and since the stability 
constraint only requires CFL< 1 this approach can potentially be quite efficient. 

Secondly, we take into account the multidimensional correctio ns required to balance th e 
V ■ B terms in a form that is simpler t han originally proposed by lGardiner fc Stond ( l2008l ). 
and analogous to lCrockett et al.l (|2005[ ). Our tests suggest that the accuracy and robustness 
of the algorithm are not affected by this simplification. 



Finally, using in p articular the id eas in iBerger fc Colellal ( 119891 ) for adaptive-mesh- 
refinement (AMR) and iBalsaral (120011 ) for the divergence-free coarse-fine interpolation of 
the magnetic field in newly refined grid patches, we have implemented an AMR version 
of our CT U + CT MHD scheme. The code in which the implementation is carried out 
is CHARM (IMiniati fc Colellal l2007al ). It includes various other physical modules, namely 
self-gravity, collision-less d ark matter particles and c osmic expansio n for cosmological a ppli- 
cations, radiat ive cooling (jMiniati fc Colellal l2007bl ) . cosmic-rays (iMiniatil l200ll . 120071 ) and 
dust particles dMiniatilljom 



The remainder of this paper is organized as follows. The algorithm is discussed in detail 
in Sec. [2J In Sec [3] we present results for an extensive set of problems that test the accuracy 
and robustness of the code. These tests include a convergence study in smooth flows, a suite 
of Riemann problems in one and two dimensions, the Orszag-Tang Vortex as well as the 
rotor problem, carried out on a uniform grid. Finally, extension to AMR and to the case of 
cosmological applications are described in Sec. S] and El respectively. We have tested these 
extensions with a problem involving the interaction of a magnetized interstellar cloud with 
a strong shock, and the formation of a galaxy cluster in a fully cosmological simulation. The 
paper closes with a brief summary in Sec. El 



2. Numerical Scheme 

In this section we first provide an overall description of the full CTU + CT algorithm 
and, following that, we discuss the implementation details. We begin by a description of the 
space discretization, data structure and notation. 
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Fig. 1. — Control volume and discrete representation of physical quantities. 

2.1. Preliminaries 

2.1.1. Discretization, Variables and Operators 

The algorithmic operations are carried out on a discrete representation of the continuous 
D-dimensional space given by the cubic lattice, i = {io, ■■■,i-D-i) G iP. The computational 
domain, referred to as a grid F, is a bound subset of iP and provides a finite-volume 
discretization of the continuous space into a collection of control volumes, faces and edges. 
Each control volume is identified by an index i = [io, ■■■,i-D-i) G F and corresponds to a 
region of space, 

Vi=[ih,{i + v)h], (5) 

where h is the mesh spacing, and v = (1, 1) is the vector whose components are all equal 
to one. The face-centered discretization of space based on the same control volumes is: 
Ff'* = { j ± ie'^ : i G F}, where e'^ is the unit vector in the d direction. Ff'* indexes the cell 
faces of F normal to e'^ representing the areas 

A^,. = [ii + e')h,ii + v)h], i + ^e'^erf. (6) 

Finally, the edge-centered discretization of space is: Ff'' = {i± ^e'^'^ ± |e'^2 ■ ^ y, d di 
^2}. F^'' indexes the edges of the cells in F aligned with e*^ representing the lengths 

^i+i(e^i+e^.) = [ih + ^(e^^ + e'^^), (i + v)h], i + ^(e'^^ + e^^) G Ff . (7) 
Fig. [1] illustrates a control volume with the various types of discretization. 
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Given the above discretization, we define cell-centered discrete variables on T as 

(f):T 

and denote by (pi G M™ the value of at cell i G F. Similarly we define face-centered vector 
fields on as 

F = (Fo,...,Fd_i) ,F,:rf ^M'", 

and denote by -F^ j+igd G M*" the value of at i + le*^ G Tf'^ , and also define edge-centered 
vector fields on r^'' as 

and denote by E^^-^i^^a^^^a^-^ G M™ the value of at i + l{e'^^ + e'^^) G rf . Finally, for 
face-centered fields we introduce the discretized divergence operator 

1 

■ = /I E (^^.i+|e^ - ^<^,i-|e^) e r, (8) 

* d=0 

and for edge-centered fields the discretized curl operator 



d\,d2 



Time is also discretized into a number of finite intervals of variable size. In particular, 
we write t^^^ = + At"', where t indicates the solution time, n indicates the integration 
step, and At" the time-step interval. 



2.1.2. Time Integration 

The coupled system of equations ([I])-(j4j), with the allowance for a non-zero source term 
S{U), can be cast in the following compact form 

dU 

^ + V.F = S, (10) 
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where the conservative variables, U, and the associated fluxes along the direction d, Fd, are 
defined as 

/ P \ / PUd \ 



U 



pUo 

pe 
Bo 



pur>-iUd + P 5dr>-i - ^D-i^d 
Ud {pe + p) - BdB ■ u 
uoBd - BoUd 



(11) 



\-Bd-i/ \ un-iBd - Bn-iUd ) 

Following Godunov's approach and its higher order extensions, we can conveniently formulate 
a numerical integration scheme based on the conservative properties of the system ffTOj) . 
In this approach one follows the evolution of cell-centered volume-averaged conservative 
variables, defined as 



f/(t", X) dV. 



;i2) 



i JVi 



The evolution equat ion is obtained upon suitable manipulation of the original continuous 
Eq. (fTOl) . and reads ( jLevequd 119981 ) 



where the face-centered time-averaged fiuxes along the direction d are defined as 

1 



(13) 



d,i+\e'>- 



dt 



Fdit.x) dA. 



(14) 



If the source term is non-stiff, we can obtain a second order estimate for it by the simple 
time-average, ^"+^ ~ K^'^ + ^"+^). 

Note that given the fiux along a direction di, -F^^^.,. igdi , the component corresponding 
to the magnetic field along the direction d2 di, defines a face centered electric field along 
d according to 

Ed^i+i^di = -SddidiF'dui+^e'ii^^^' (-'-^) 

where we use subscripts to indicate directions and centering, and superscripts for compo- 
nents. Indeed, Godunov's scheme also updates in time the cell-centered magnetic fields 
variables. As already pointed out, however, the updated magnetic field in general does not 
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remain solenoidal. Therefore, we adopt instead a CT discretization strategy in which the pri- 
mary description of the magnetic field we will be using face-centered area-averaged variables 
defined as 

(16) 



2 ^i+le'i J A. 



Note that an estimate of the cell-centered magnetic field variables is still needed in order to 
construct the fluxes in (I14p . We will return to this point shortly. The evolution equation for 
the face-centered magnetic field variables is obtained again from a manipulation of Faraday's 
law and reads 



B 



n+l 



E' 



-{E 

+S' 



2 - ' 2 

2 fp 2 



(17) 



with (i 7^ (ii 7^ ^2,0 < d,di,d2 < 
formally defined as 



D. The edge-centered time-averaged electric field is 



p"+2 
d,i+^edi+^ed2 



1 



At 



i+e'*2) Jt 



m + l 



dt 



Ed{t,x)dL. 



(18) 



e<*2) 



and, as above, the time-centered estimate of the non-stiff source term for the face-centered 
magnetic field components is obtained by simple arithmetic averaging. The cell-centered 
magnetic field variables are then defined in terms of the pri mary face-centered values using 



the following second-or der accurate reconstruction scheme (iRyu et al.l Il998l : iBalsaral 12001 
Gardiner fc Stonelboosh 



B, 



dA 



dM 



(19) 



An important part of CT schemes concerns the calculation of the time-averaged, edge- 
centered electric fields e ntering the time update ( fT7|) . In general one employs some type 
of bilinear interpolation. iDai fc WoodwardI (jl998[ ) interpolate in space and time the cell- 
centered velocity and magnetic field at time with the solution from the Godu nov scheme 
at tina e to obtai n time- and corner-centered electric fields. On the other hand, iRyu et al. 
( 119981 ) and IBalsaral ( l200ll ) take advantage of Eq. ( 1TB]) and interpolate the face-centered 
variables th at allow reconstru ction of electric fields at cell edges. The interpolation scheme 
proposed in iRyu et al.l ( 119981 ) has the property that for plane-parallel configurations their 
multidimensional scheme reduces to th e one-dimensional schern e. This same property is 
shared by the upwind scheme proposed in lGardiner &: Stond (120051 ) which is adopted here and 
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is described in more detail in Sec. 12. 2.61 This property is important because it guarantees self- 
consistency between (a) the electric fields used for the time-update of face-centered magnetic 
variables, (b) the MHD fiuxes used for the time-update of the cell-cent ered variables and 



c) th e synchronization step of the magnetic variables given in Eq. f[T^ (iGardiner &: Stone 



20051). 



Note that by applying the divergence operator in Eq. ([8]) to the face-centered magnetic 
field evolved according to Eq. ( 1T7|) , one finds that the divergence of the field does not change 
in time. So the magnetic field remains solenoidal if initially so. This is the property of the 
CT scheme. 

The accuracy and stability of the numerical solution depend principally on how the time- 
averaged fiuxes and electric fields entering Eq. f[T^ and f lT7|) . respectively, are computed. In 
the following sections we shall describe the algorithmic details characterizing such calculation. 



2.2. Algorithm 

The algorithm described in this section computes second-order accurate face centered 
fiuxes and edge-centered electric fields required for the update of the cell-centered fiuid 
variables and face-centered magnetic f ield variables in Eq.(IT^ and (fT7|) . respectively. It uses 
a combination of the CTU algorithm (jColellalll990t ISaltzmanlll994j ). and the CT scheme, as 



m 



Gardiner fc Stond ( l2008l ). Unlike these authors, though, for the reasons indicated in the 



Introduction we will employ the full CTU scheme with 12-solve per cells (in 3-dimensions). 

Assuming the solution at time t", f/" to be known, the outline of the algorithm is as 
follows: 



1. transform from conservative to primitive variables, f/f ^ PFj", and synchronize cell- 
centered and face-centered magnetic field values 

2. compute limited slopes along the coordinate directions, SW^^^. 

3. do characteristic tracing to extrapolate in space and time primitive variables from cell 
centers to cell faces, Wi^±^d- Also include effects of source term here. 

4. convert back to conservative variables, PFi,±,d ^ Ui^±^d- 

5. apply corrections to Ui^±^d due to transverse gradients and obtain multidimensionally 
correct time-averaged fluxes F'^\^ , and electric flelds E^^^ ^ ^ 



- 11 - 



6. update primary variables, f/" ^ f^""*"^ ^ind -B^^i^d <— B'^^l^^ and synchronize cell- 
centered and face-centered magnetic field values. 

7. update solution time t ^ t + At and compute new timestep according to the Courant- 
Friedrichs-Lewy (CFL) condition for stability, 

where the CFL number is, Ccfl < 1, is the d— component of the velocity field, c/ is 
the fast magnetosonic wave speed and the max value is computed over all directions d 
and over all cells in the computational domain. 

In the following subsections we provide additional details about the algorithm. 



2.2.1. Primitive Variables 

The first part of the algorithm consists in reconstructing the numerical solution within 
the spatial domain of the mesh. This requires estimating gradients and eventually extrapo- 
lating state variables in time and space from cell centers to cell faces with the use of char- 
acteristic tracing. Such operations are usually done in primitive space, where the variables 
are defined as 

W = {p,uo, . . . ,ur,-i,pg, Bq, . . . , Bn^i)'^ , (21) 

as it simplifies the characteristic analysis of the system. The evolution of the system in terms 
of primitive variables is given by a system of equations in non-conservative form, 

where = VuW ■ S and Ad = VuW ■ VuFd ■ VwU, and VuW and VuW are the Jacobian 
of the transformation from primitive to conserved variables and vice-versa. Given their 
importance for the construction of a sound numerical sche me, the properties of the operators 



Ad h ave been studied in great details in the literature (e.g. iBrio fc Wulll988l : iRoe fc Balsara 



19961 ). Here it suffices to make the following observations. In one dimension, d = 0, the 
parallel component of the magnetic field is constant, Bq = const., and Aq effectively becomes 
a 7 X 7 matrix. In general the operator Aq is characterized by 7 eigenvalues, and associated 
left and right eigenvectors, with values u, u ± Cs, u ± ca and u ± Cf obeying the hierarchy: 
Cs < Cyi < c/. The first eigenvalue listed above corresponds to the usual entropy wave, and the 
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other six to the three pairs of MHD waves (slow magnetosonic, Alfven and fast magnetosonic) 
propagating downstream (+) or upstream (-) in the flow, respectively. Because up to 5 of 
the eigenvalues may actually coincide, the system is not strictly hyperb olic, so care must 



be taken to avoid siri gularities with expressions involving the eigenvectors (IBrio &: Wulll988 



Roe fc Balsaralll996l ). Finally, in more then one dimensions, because Bd is not affected by 



gradients in the d direction, it has been customary to use the one- dimensional analysis when 
formulating the predictor step for higher-order Godunov-like MHD algorithms. However, 
that leads to neglect of terms oc dB^/dxa, which are not necessarily null in multidimensions. 
In fact, it is important to include these te r ms to avoid degrading the solution accuracy, as 
recently pointed out by iGardiner fc Stond (120051 ). 



2.2.2. Slopes 

After the conversion from conservative to primitive variables, ^ ^^i^i ) 
synchronization of cell-centered and face-centered magnetic fields according to Eq. (1191) . we 
proceed to the calculation of the slopes along each direction d as follows. First, central and 
side slopes are estimated, respectively, as 

6,,oW, =^ {W-^^, - Wl^,) , (23) 

Sd,-Wi=W,^-Wl^„ (24) 
6,,^Wi =Wl^^, - W^r, (25) 

and then limiting is applied component-wise either in primitive or in characteristic space. 
We use van Leer's limiter defined as 



5" ((5o,5_,5+) 



sgn((5o) min(|5o|,2|(5_|,2|5+|) if > 
otherwise. 



In the case of primitive limiting the limiter is applied directly to each component k of the 
primitive slopes (125]) i.e. 

= S^\S,,oW,', 6d,-W,', 5rf,+iy^). (26) 

In the case of characteristic limiting the limiter is applied to the components of the primitive 
slopes in characteristic space, namely 

SdWi = ^«'r^ (27) 

k 

= 5''^(ag,a^,a^), (28) 
a# = I'-Sd^^Wi, # = 0,+,-, (29) 
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where, = /^(Ty/^) and r'^ = r^lWl^) are the left and right eigenvectors of the operator Ai^a 



2.2.3. Normal Predictor 



Next we extrapolate the primitive variables from cell centers to face centers along each 
coordinate direction, d, by taking into account the time-averaged effect due to the slopes 
just computed. This is done most conveniently by using the 1-D versions of the MHD 
equations in primitive form. There is no evolution in the (i— component of the magnetic 
field due to derivatives along the direction, so the normal components of the magnetic 
field on cell faces are straightforwardly provided by the face-centered component of the 
magnetic field, without need for even g eometrical extrapolation. On the other hand, as 
pointed out by iGardiner fc Stond ( l2005l ) in multidimensional MHD the reconstruction step 
must include terms proportional to dBd/dxd terms (no summation over repeated indices is 
implied here). These terms arise from the requirement to balance the divergence terms in 
more than one dimension and their neglect can cause serious degradation of the numerical 
solution. The reconstruction of the primitive variables onto cell faces can be done with 
various degrees of accuracy. Typically, one uses a piecewise-linea r (PLM) or piecewise- 
parabolic (PPM) reconstruction scheme ( IColella fc Woodwardlll984j ). Although we mostly 
use a PPM algorithm for the tests presented here, for simplicity we illustrate the case of 
PLM reconstruction. So in this case the extrapolation of the primitive variables in space 
and time from cell centers to faces along the direction d takes the form 



B 



d,i,±,d J 



^dA+e'i 



f±I 
I 




P±{5dWl 



+ 



At S. 



■n,MHD^ 
d,i 



(30) 



where we have explicitly separated out the reconstruction of the normal component of the 
magnetic field (the hat symbol in the notation indicates that the components corresponding 
to the normal magnetic field are omitted). In addition, we have used the projector operator 
defined as 

P±{W) = 5^ (4 ■ W)rk (31) 

where are eigenvalues of Ai^a- This projector operator filters out the components of the 
gradients that propagate away from the cell interface. However, when a Riemann solver of the 
HLL family is employed, in order to obtain second-order accuracy the filter is switched off and 
both in the PPM and PLM cases, and the summation is carried over all waves, irrespective of 
their sign (this is further discussed in Sec J2.2.5|) . Finally, S'd,MHD represents the MHD source 
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term required in multidimensional MHD which we implement in the form ( ICrockett et al. 
2005h 

/ 

Bq 
P 



S. 



n,MHD 
d,i 



\ 



Bo- 



P 

B u 



dxd 



(32) 



where, as usual, d ^ d\ ^ d2 and, in this particular case, < cii < < D- The normal 
predictor step is completed by the final corrections for a non-stiff source term according to 

At 



(33) 



2.2.4- dT Extended Corner Transport Upwind 

After the primitive variables have been extrapolated to cell faces, we add corrections 
due to gradients parallel to the cell faces. We find it convenient at this point to convert 
back to a conservative representation, thus we operate the transformation: W7^_|_ ^ ^ U^^ ^. 
Following the CTU scheme the corrections are expressed in terms of transverse flux gradients, 
with fluxes obtained from a Riemann solver, TZ. In accord with the CT scheme, however, 
for the update of magnetic field variables we use gradients of edge-centered electric fields 
suitably interpolated in space and time from their face and cell centered values. Both the 
interpolation procedure, X^, and the Riemann solver, TZ, will be specified at the end of this 
section. It suffices here to say that the time centering and interpolation accuracy of the 
interpolated edge-centered electric fields is consistent with that of the MHD fiuxes returned 
by the Riemann solver. 

The steps involved in the modified CTU update can then be summarized as follows: 

1. Use a Riemann solver, TZ{U^^^^ , jjRight^^ -(.q obtain a first estimate of the fiuxes across 
cell faces along each direction, d 

FS^.'^ = mi,+4^Ul^^,-,d,d). (34) 

2. Use the newly obtained fiuxes with the primitive solution at time to interpolate the 
electric fields from cell-faces and cell-centers to cell-edges 

-E'i°+iedi + igd2 = 2:£;(F]°^^i^di , +ied2 , W^), (35) 
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where the * symbol indicates that the interpolation requires values of the arguments 
at various cell centers and faces. 

3. As part of the CTU prescription to obtain (1, 1, 1) diagonal coupling, apply corrections 
to density, momentum and energy components of Ui^±,di, due to one set of transverse 
flux derivatives along d2. For each face di, there will be D — 1 such corrected states, 
one for each direction perpendicular to di. The corrected states read 

At 



4. Likewise, use the CT scheme to correct magnetic field components affected by the same 
set of transverse flux derivatives. There are D — 1 magnetic field components of U^_^ 
on the face di which are affected by the transverse flux along d2- First, the component 
along di , which we indicate with Ub^^ , is corrected as 

TT^ — TT"^ — c ( _ plD 

'^Ba^,i±,dud2 -^Bdi,i±,di °'i'^3d22^^ y-^d3,i±^ed'i + ^e'i2 ^ds,i±^e'ii-^e'i2 

Second, we correct the magnetic field component, Ub^,^, p arallel to the cell face an d 



directed along the remaining direction d^ ^ d2 ^ di. As in iGardiner &: Stond (120081) 



for this component we average the contributions from the faces at i + \e'^^ and i — ^e 
obtaining 

At 



TTn _ Tjn _ 

Ba^,i±,di,d2 - ^Sd3,i±,di '^i'^3d2g^^ 



I / T-ilD 7-.1D 
+ I -^di,i-ie''3 + ied2 ^di,i-^e'i3-.^e''2 



{3t 



5. Use a Riemann solver to obtain fluxes for each pair of states corrected for transverse 
fluxes. This provides D — 1 fluxes per cell face. 



dl ^ d2, < di, 6/2 < D 

6. Obtain new interpolated values of the electric field from the averages of the above 
computed fluxes 

-Erfi+igdi + igda = ,,+ igdi,F^2,*+|erf2, (40) 

where 

^dl,i+ie'*l = ■^{F^^,i+le'il,d2 + ^di,i+^e'il,d3) (41) 
1 

-^d2,i+fe'*2 = 2(-^rf2,i+fe''2,di + -^d2,i+|e''2,d3)- (42) 
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7. Compute final corrections to the density, momentum and energy components of Ui^±^d 
due to transverse fluxes using the above Riemann solutions, according to 



At 

2Ax 
At 



2Ax 

d ^ di ^ d2, < d, di, < D 



(43) 



Likewise, use the CT scheme to apply final corrections to magnetic field components. 
In analogy to Eq. fl57|) in step 4 we write 



Bi,i,±,d — '-^Bd,i,±,d ^ddid2 



At 

2A^ 



E. 



di,i±\e<l+};e'^2 



E. 



di,i±\e<i-\e'^2 



(44) 



Finally, in analogy to Eq. (13811 the magnetic field components parallel to the cell face 
are updated according to the CT scheme as 



TT ~'~2 — TT^ c 

'~^Ba^,i±,d — '-^Ba^,i±,d ^ddid2 



At 

2Ax 



E 



+ E, 



^d,i-ie''l + ie''2 



~E. 



d,i-^ed-i-^e''2 



- E 



- E 



'^d2,i-^e'^i- + le^ ^■ 



d2,i-'^e'^i-^e^ 



(45) 



9. Compute final second-order estimate of time averaged fluxes 



i+e''',-,d^ 



d). 



(46) 



10. Compute final estimate of the electric field using the above fluxes and an updated 
value for the cell-centered conservative variables using the averaged fluxes in Eq. ( HTl) . 
namely 

(47) 



771"'+ 2 7- /'p"'"*"2 fr"~'~2 \Al^\ 



where 



D-l 



WuW ■ J2 {Pd 



d=0 



and the operator symbolizes the transformation from conservative to primitive 

variables. 
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11. Finally update cell-centered conservative variables using Eq. (1131) with the fluxes in 
fH6|) and the face-centered magnetic field variables using Eq. f|T71) with the electric 
field in fH7|) and synchronize the magnetic variables using Eq. flTI?]) . 



2.2.5. Riemann Solver 

For the purpose o f this p aper, we have implemented the HLLD solver recently developed 
by lMiyoshi &: Kusand ( 20051 ). This is an extended version of the original HLL solver for the 
MHD, which includes the entropy, Alfven and fast magnetosonic waves. The solver appears 
quite accurate and robust and relatively inexpensive. As already pointed out in order for the 
fluxes returned by HLL-type solvers to be second-order accurate, the projector P is modified 
in such a way that the summation is carried over all waves, irrespective of their sign. This is 
because the fluxes returned by these solvers are built on the spatially reconstructed solutions 
on the left and right hand sides of the cell interfaces. The solver is extensively documented 
in the original paper and its description will not be repeated here. 



2.2.6. Interpolation Scheme for the Edge-Centered Electric Field 



In this section we describe the scheme used to interpolate the face-centered electric 
fields returned by the Riemann solver onto cell edges. In Sec. 12.2.41 this was indicated with 
the notation, Xe- It is important for the stability of the overall algorithm to choose this 
interpolation scheme in such a way that there is consistency between the time-update of 
the face-centered magnetic field, the cell-centered variables and the synchronization step of 
the magneti c variables. Eg. (flQl). For this purpose we have adopted the upwind scheme 
described in iGardiner &: Stond ( l2005l . l2008l ). This will be described next, for the case of 
three dimensions. In this case each cell edge is shared by four adjacent faces from which the 
electric field can be interpolated. Hence, these four interpolated values will be arithmetically 
averaged. The scheme for Te{) uses three arguments. 



E. 



A2 , 



where the * symbol indicates that the interpolation requires values of the arguments at 
different cells and faces. The face-centered fluxes define the face-centered electric fields 
according to Eq. (ITSl) . and the cell-centered primitive state variable, Wi, is used to define a 
cell-centered electric field according to the usual formula 



E. 



^ddid2 
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with the velocity and magnetic field variables given by the components of Wi. The cell- 
centered electric field is used together with face-centered electric fields f|T5|) to define the 
following transverse quasi-cell-centered gradients 



2— 



E. 



dA 



'^1 / i+4e'*i 



(4J 



In turn, these gradients are used to define quas i-face-centered gradients of the electric field 
by the upwind scheme (iGardiner fc Stonell2005l ) 



5E, 



5x, 



d2 



'^2 / i+ie''2 



SEg 

5x, 



''2 / i+fe'^a 



+ 



SEg 

5x, 



2 / i+|e<*2 



if Ui^i^d2 > 



if Mi+le''2 < 



otherwise. 



With the above definitions, the interpolated edge-centered electric field is defined as 



E 



d,i+§e'^i + |e'*2 



+ 



+ 



1 
4 

Ax 



Ax 



^d,i+iedl + E^ 



d,i+e''2 + |e 



i+|e''i+|-e'*2 



SE, 

SXd2 



.1 

SE, 
6xd2 

SEd 

5x, 



d,i+ied2 + E, 



d,i+e''-i + \e<^2 



i+|e'*i+|e''2 



(49) 



3. Tests 

3.1. Convergence Rates in Smooth Flows 

In this section we test the correctness of our implementation by measuring the conver- 
gence rate of the numerical solution returned by the code. The test is based on the propa- 
gation of Alfven, fast and slow MHD waves. The waves have small amplitude, 5 = 10~^, so 
that we expect to observe the nominal second order convergence rate predicted by numerical 
analysis. 

In the following tests the initial conditions are provided for the primitive variables by 
defining the unperturbed state, W, and the superposed perturbation, 6W, corresponding 
to the wave. The size of the computational box is, L = 1, the geometry is one- or two- 
dimensional and the boundary conditions are periodic. The adiabatic index is 7 = 5/3. 
While we have experimented with different choices of orientation of the wave-vector with 
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respect to the grid, namely k/27r = (1,0), (0,1), (1/^2,1/72), (2/75,1/^5), we find 
the same convergence rates in all cases. Thus, we simply report the results for the few 
representative tests listed in Table [H 

In order to measure the rate at which the numerical solution converges, for each problem 
we carry out a set of 5 simulation runs employing N^eii = 16, 32, 64, 128, 256 for a total range 
of 16. For each run the time-step size is fixed, scales inversely with Nceiu and corresponds 
roughly to a CFL number 0.8. The convergence rate is measured using Richardson extrapo- 
lation. Given the numerical solution at resolution r we first estimate the error at a given 
point («, j), as 

er;i,j = gr(^, j) " gr+l(^, j)> (50) 

where g^+i is the solution at the next finer resolution, spatially averaged onto the coarser 
grid. We then take the n-norm of the error 

-^n ~ ll^rlln — {^,\^'r;i,j\ '^ijj ) (51) 

where, Vij = Ax^ is the cell volume, and estimate the convergence rate as 

_ ln[L„(£,)/L„.(£,)] 
"~ ln(Ax,/Ax,) ■ ^ ' 

For each case listed in Tabled! we produce a corresponding Table [2lll] reporting the Li, L2 
and Loo norms of the error and the corresponding convergence rates, Ri, R2 and Roo, as 
defined above. 



Table 1: Run Set 



run 


6 


k/27r 


Type 


A 


10"^ 


(1,0) 


Alfven 


B 


10-5 


(2/v^,l/v^) 


Alfven 


C 


10-5 


(1,0) 


Fast 


D 


10-5 


(2/75, 1/V5) 


Fast 


E 


10-5 


(1,0) 


Slow 


F 


10-5 


(2/v^, 1/V5) 


Slow 
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3.1.1. Alfven Waves 



The test based on the propagation of the Al fven wave is characte rized by the following 
unperturbed state and superposed perturbation (ICrockett et al.ll2005l ) 



W 



Po 



Po 

B0/V2 
B0/V2 





6W 



( \ 




-ca 






b sin(k ■ x) 



(53) 



where po = -Po = -Bo = 1 and Ux = Uy = Uz = 0, ca = B^j .^fp^ = 1, is the Alfven speed, 
and k and x are the wavevector and position vector respectively. The convergence rates for 
the perturbed quantities are summarized in Table |2l In both tests, with different k, both 
the z— components of the velocity and magnetic field converge with second-order accuracy. 
The error on the unperturbed variables (not reported) is much smaller and converges at the 
same rate as the perturbed variables until dominated by the machine round-off error. As 
already pointed out, tests with different orientation of k show the same convergence rates. 



3.1.2. Fast and Slow Magnetosonic Waves 



For fast an d slow magnetosonic waves the unperturbed state and the superposed per- 
turbations read (ICrockett et al.l 120051 ) 



W 



/po \ 

Uy 

Uz 

Po 
Bobx 

B(,by 

V / 



5W 



( 



Po 



iV2clby - 

-iV2clbx 



Poc] 



-V2B,{cl-cl)ky/cl 
V2Bo{ci - cl)kxlc\ 



\ 



I 



S sin(k 



(54) 



where, b is the unit vector along the magnetic field vector and is at 7r/4 radians with respect 
to the unit vector = k/fc, Cg = -^Z 7P0 / Po is the gas sound speed, is the speed of fast 
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or slow MHD waves and all other symbols take the meaning and values as in the previous 
section. 

The convergence rates are reported in Table [3] and H] for fast and slow waves, respectively. 
As with Alfven waves, the errors in the perturbed variables converge with second order 
accuracy, while the error in the unperturbed variables are much smaller and converge at the 
same rate until they are affected by machine precision. 



3.2. Riemann Problem 



W[xo,t = 0] 



(55) 



We next consider a set of Riemann problems from the literature which are standard 
tests for MHD algorithms. The Riemann problem is described in general by the following 
initial-value problem 

Wieft if xo < 0.5 
Wright if xo > 0.5 

where Wieft/Hght represents the primitive variables to the left/right of the initial discontinu- 
ity. The set of problems and corresponding initial conditions are summarized in Table |5l 
except for the value of the x— component of the magnetic field which will be specified for 
each problem explicitly. In all cases we use a CFL number Ccfl = 0.8, third order PPM re- 
construction scheme and characteristic limiting. The domain size is L = 1 and the boundary 
conditions are simply, W{xq = 0,t) = Wieft, and, W{xo = l,t) = Wright- 



We should point out that iMiyoshi fc Kusand ( l2005l ) have carried out extensive tests of 
their HLLD solver, which we employ in our code, particularly to compare its performance 
with that of Roe and other solvers of the HLLC family. We sha ll repeat some of those tests 
here. Note that while in most cases iMiyoshi &: Kusand (|2005[ ) used a first order accurate 
piece-wise constant reconstruction scheme, we have used the third order accurate PPM 
method, so the solution profiles in our plots appear sharper. The tests which we will perform 
below involve the full set of MHD waves, including those associated with t he slow mode which 
i s not included explicitly in the HLLD solver. However, in accord with IMiyoshi fc Kusano 
(120051), we find that in these tests the full MHD structure of the solution, including features 
associated with the slow mode, is correctly reproduced. 



3.2.1. Bno & Wu 



We begin with the Riemann problem presented in lBrio &: Wd ( 119881 ). listed at the top of 
Tabled The x— component of the magnetic field is = 0.75 and the adiabatic index 7 = 2 
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Fig. 2. — Brio & Wu shock tube problem: solution for t = 0.1 solved on a grid using 800 
zones. See Table |5] for the initial conditions. From left to right and top to bottom shown 
are, respectively: density, pressure, velocity components along the x and y axis, magnetic 
field component along the y—axis and temperature. 
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as in the original paper. Following common practice, the problem is solved on a grid with 
800 points along the x— axis and the solution is computed until time t = 0.1. The results 
are shown in Fig. |2J All the solution features are well reproduced, including, from left to 
right, a fast rarefaction followed by a slow compound wave, moving t o the left, arid a co ntact 
discontinuity, slow shock and fast rarefaction moving to the right ( iBrio fc Wdll988l ). As 
usual with shock-capturing methods, shocks are resolved with a couple of zones throughout 
the duration of the calculation, while the contact discontinuities, captured here within a few 
zones, tend to spread out over time. There are some oscillations in the x— component of the 
velocity field. These are not present if we adopt a PLM reconstruction scheme, but worsen 
if we switch from a characteristic to a primitive limiting scheme. 



3.2.2. Dm & Woodward 



The second Riemann problem listed in table is taken from iDai fc WoodwardI ( 1l998l ). 
In this case the x— component of the magnetic field is = A/\/A^ and the adiabatic index 
7 = 5/3. The problem is solved on a grid with 512 points along the x— axis and the solution 
is computed until time t = 0.2. The results are shown in Fig. |3l The initial conditions for 
this problem expose the full eigenstructure of the MHD system, as they produce three pairs 
of MHD waves traveling in opposite directions with respect to the initial discontinuity, in 
addition to the contact wave. The waves include fast and slow shocks responsible, among 
others, for the jumps in the pressure and velocity fields, the contact wave which appears 
in the density field alone, and the rotational discontinuity which affects the magnetic field 
components alone. As in the previous case, while all discontinuities are well reproduced, 
shocks are the sharpest features captured with about two zones. 



3.2.3. Ryu & Jones 



The third Riemann problem in our table[5]is one of the many problems studied by lRyu fc Jones 



(119951 ). Here = 0.7 and 7 = 5/3. The problem is solved on a grid with 512 points along 
the X— axis and the solution is computed until time t = 0.16, as in the original paper. The 
results are shown in Fig. HI The solution to this problem includes, from left to right, a 
hydrodynamic rarefaction, a switch-on slow shock, a contact discontinuity, a slow shock, a 
rotational discontinuity and a fast rarefaction. Switch-on/off fast /slow shocks are interesting 
MHD solutions that cause the tangential magnetic field to turn on/off, respectively. As in 
the previous tests the structure of the solution is well reproduced, including the features 
associated with the slow MHD mode, and the discontinuities are captured within a few to 
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Fig. 3. — Dai & Woodward shock tube problem: solution for t = 0.2 solved on a grid using 
512 zones. See Table E] for the initial conditions. From left to right and top to bottom, 
shown are, respectively: density, pressure, energy, three velocity components magnetic field 
components, along the y and z axis, and e = tan-^{B,/By). 
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Fig. 4. — Ryu & Jones shock tube problem: solution for t = 0.16 solved on a grid using 
512 zones. See Table E] for the initial conditions. From left to right and top to bottom, 
shown are, respectively: density, pressure, energy, three velocity components, magnetic field 
components along the y and z axis, and e = t&n-^{B,/By). 
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several zones. 



3.2.4- Fast rarefaction 



The last one-dimensional Riemann problem listed in table [5] is similar to the ones pre- 
sented in lEinfeldt et al.l ( 1l99ll ). Its purpose is to test the method/code robustness in the case 
of flows in which the energy is dominated by the kinetic component and unphysical states 
with negative density or internal energy can arise. This problem is solved using a grid with 
256 grid points. The solution at time t = 0.1 is shown in Fig. |3 The results in Fig. Oshow 
a stable solution which correctly reproduces the two rarefaction wave s propagating awa y 
from the grid midpoint. A very similar test has al so been performed by lStone et al.l ( l2008l ). 
with similar results, and lMiyoshi fc Kusand (120051 ). The latter authors actually use a faster 
expansion velocity, namely l = =p3, in their initial conditions and show that when used 
with a (first order) Godunov's method the numerical solution remains stable. While we are 
able to reproduce their result we find that at high resolution some spurious oscillations can 
be generated when a higher-order Godunov's method is used. This may suggest the need for 
artificial viscosity in the case of highly supersonic expansions with higher-order Godunov's 
methods. 



3.2.5. Inclined Dai & Woodward Shock-Tube Problem 

We have repeated the problem in Sec. 13.2.21 with the initial discontinuity inclined with 
respect to the grid such that its normal has components fi = (2, l)/\/5. This problem tests 
the ability of the code to reproduce one-dimensional solutions when they are not aligned 
with the grid. Besides the numerical tests, there are a few other complications related with 
the numerical set-up of this problem. First, the boundary con ditions need revision. To avoid 



the implementation of "shift-periodic" boundary conditions (ITothl 120001 ) . we use a domain 
with size {2L, L)2/\/5. This allows us to accommodate two adjacent identical problems 
along the direction, n = (2, l)/^/5, and apply periodic boundary conditions to the domain. 
In addition, since the magnetic field rotates across the discontinuity, the initialization is 
non-trivial and, unless a special precaution is taken, e.g. by deriving the magnetic field 



^Obviously, since the sequence of states will be Wie/t, Wright, Wieft, Wright, the interaction at the second 
interface will produce a similar Riemann problem with inverted initial conditions. However, we will present 
only part of the solution, selected for proper comparison with the analogous one-dimensional problem in 
Sec. [3T2l 
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Fig. 5. — Fast rarefaction: solution for t = 0.16 solved on a grid using 512 zones. See Table [5] 
for the initial conditions. From left to right and top to bottom, shown are, respectively: 
density, gas pressure, x— component of the velocity and y— component of the magnetic field. 
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Fig. 6. — Inclined version of Dai & Woodward shock tube problem: solution for t = 0.2 solved 
on a grid using 512 zones. See Table [5] for the initial conditions. From left to right and top 
to bottom, shown are, respectively: density, pressure, energy, three velocity components and 
three magnetic field components. 
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from a vector potential, the initial magnetic field will not be divergence-free (this is not 
an issue in the special case in which n is along the diagonal and there is symmetry among 
the coordinate axis). In our case no such precaution is taken and we just remap the initial 
conditions onto the rotated grid. As a result there is a jump in the normal component of 
the magnetic field across the discontinuity. This causes some minor artifacts with respect 
to the one-dimensional solution. This is acceptable since we are interested in making sure 
that the structure of the one-dimensional solution is reproduced with fidelity and the waves 
propagate at the correct speed. 

In order to solve the problem, we have covered the domain with a grid of 1144 x 572 cells. 
This corresponds to 512 cells along the direction perpendicular to n, which is equivalent to 
the resolution used in Sec. 13.2.21 The results are shown in Fig. E] where we plot the values 
of the solution along the first row of the computational domain, starting from |L + 6l to 
+ The starting point is shifted by 6l ~ 0.2 to the left, to hide the perturbation of Wiejt 
due to the interaction with the state to its left (which is Wright)- The vectorial components 
(m, V, w) are the equivalent of (x, y, z) in the rotated system. 

Fig. [5] shows that the one-dimensional solution is correctly recovered when the plane of 
the discontinuity is inclined with respect to the grid. We notice some osc illation at the fast 



magne tosonic shocks, which to some extent are also present in Fig. 7 of iGardiner &: Stone 



(120081 ). These are probably associated with the oscillations in the normal component of the 



magnetic field, also reported in the last panel of Fig. El These can perhaps be suppressed 



with more aggressive limiters ( iLondrillo fc del Zannall2004l ). There is also a spurious feature 



ahead of the left-moving rotational discontinuity, which is probably due to the non-solenoidal 
character of the initial magnetic field. However, none of these features affect either the jump 
conditions or the wave speeds, as attested to by the very good correspondence between the 
solution in Fig. E] and the one- dimensional counterpart in Fig. 121 



3.3. Magnetic Loop Advection 



In this section we test the ability of the code to follow the advection of a loop of weak 
magnetic fie ld frozen in a background flo w. This problem is non-trivial for conservative 
schemes and lGardiner fc Stond (l2005l . l2008l ) emphasize that spurious results will be produced 
as a result of improper account of the MHD source terms entering the predi c tor st ep, as 
discussed in Sec. 12.2.31 The initial conditions are detailed in iGardiner fc Stone J2008I ). The 
loop is basically a tube of magnetic flux frozen in a medium with unit density and pressure, 
P = Pgas = 1 and uniform advection velocity, uioop, to be specified below. 
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We carry out the test both in a 2D or a 3D geometry. For the 2D case, we ahgn the 
loop axis with the X2 coordinate axis of the computational domain. The vector potential 
from which the magnetic field is initialized is then given by, A = (0, 0, A2), with 

_ r Bo{R-r) if r<R, 
"^'-[o if r>R, 

where Bq = r = a/xq + x^, and, R = 0.3, is the radius of the tube. The computational 
domain itself consists of a rectangular box of dimensions (2A^, A^) with periodic boundaries, 
and the loop advection velocity is uioop = (2, 1). 

For the 3D configuration, the axis of the tube is tilted around the xi axis by ^ = arctan 2 
radians (clockwise) and the tube is advected with velocity uioop = (1,1,2). The domain is 
still periodic but has dimensions [N, N,2N). The vector potential is still defined as in fl56l) 
but with respect to the new coordinates 

x'q = xq cos 9 + X2 sin 9, 

x'l = xi, (57) 
x'2 = —Xq sin 9 + X2 cos 9. 

The results of this test are first illustrated in Figures [7] and M which show the magnetic 
pressure at t = 2 for the 2D and 3D cases respectively. In the former case = 64, while 
in the latt er case N = 128. The re s ults a re comparable to other MHD implementations, in 



particular iGardiner fc Stond (l2005l . |2008| ). As pointed out by those authors, the magnetic 
pressure suffers dissipation mostly close to the loop center (where the curl of B is singular) 
and boundary. However, both in the 2D and 3D cases, the loop retains to a good extent its 
initial symmetry and energy. This latter point is further illustrated in Fig. [9l reporting the 
evolution of the normalized magnetic pressure up to t = 2 for three different resolutions. One 
important aspect of the 3D version of this test is to check the ability of the scheme to keep 
the magnetic field component along the loop axis, -B3, close to zero. This is importan t here 



because we have not strictly followed the recommendation of iGardiner fc Stond (120081 ) when 
implementing the MHD source terms in the predictor step. The evolution of S3 is reported 
in Figure [TOl again for three different values of the resolution, up to t = 2. Due to the simpler 



form of the employed MHD source terms, (I-B3I) is larger than found in IGardiner fc Stone 



(I2OO8I ) at the same time (t = 1), but only slightly so and S3 remains negligible compared to 



the total magnetic field. 
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Fig. 7. — Loop advection in 2D: magnetic pressure at t — 2 from a calculation using A'^ = 
and corresponding colorbar (top left). 
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Tal)le 2: Convergence Rates for Alfven Waves. 



Ncells 



Ri 



R2 



Lr 



Rr, 



L 



R^ 



Ri 



Rr, 



Case K:5 = 10 



16 

32 
64 
128 



z— velocity 
7.5E-09 - 

1.9E-09 2.0 
4.8E-10 2.0 
1.2E-10 2.0 



1.7E-08 

4.2E-09 
l.lE-09 
2.7E-10 



2.0 
2.0 
2.0 



4.7E-08 
1.2E-08 
3.0E-09 
7.5E-10 



2.0 
2.0 
2.0 



-\ k = 2^(1,0) 
z— magnetic 
7.5E-09 - 
1.9E-09 2.0 
4.8E-10 2.0 
1.2E-10 2.0 



1.7E-08 

4.2E-09 
l.lE-09 
2.7E-10 



2.0 
2.0 
2.0 



4.7E-08 

1.2E-08 
3.0E-09 
7.5E-10 



2.0 
2.0 
2.0 



Case B: (5 = 10' 



16 
32 
64 
128 



z— velocity 

7.6E-08 - 
2.0E-08 1.9 
4.9E-09 2.0 
1.2E-09 2.0 



1.2E-07 
3.1E-08 
7.8E-09 
1.9E-09 



2.0 
2.0 
2.0 



2.4E-07 
6.2E-08 
1.5E-08 
4.0E-09 



-\ k = 27r(2, l)/v^ 
z— magnetic 



1.9 
2.0 
2.0 



7.8E-08 
2.0E-08 
5.0E-09 
1.2E-09 



2.0 
2.0 
2.0 



1.2E-07 
3.1E-08 
7.8E-09 
2.0E-09 



2.0 
2.0 
2.0 



2.3E-07 
6.2E-08 
1.6E-08 
4.0E-09 



1.9 
2.0 
2.0 





Fig. 8. — Loop advection in 3D: cut of the magnetic pressure at 2; = 0.5 of the magnetic 
pressure t — 2 from a calculation using = 128, and corresponding colorbar (top left). 



-33- 



Table 3: Convergence Rates for Fast MHD Waves. 



Ncells Li Ri L2 R? ioc -Roc Li Ri L2 R2 -too Rc 



Case C:6 = 10-^ k = 27r(l,0) 
density-gas x— vel-gas 

16 7.5E-09 - 1.7E-08 - 4.6E-08 - l.lE-08 - 2.4E-08 - 6.6E-08 

32 1.9E-09 2.0 4.2E-09 2.0 1.2E-08 2.0 2.7E-09 2.0 6.0E-09 2.0 1.7E-08 2.0 

64 4.8E-10 2.0 l.lE-09 2.0 3.0E-09 2.0 6.8E-10 2.0 1.5E-09 2.0 4.3E-09 2.0 

128 1.2E-10 2.0 2.7E-10 2.0 7.5E-10 2.0 1.7E-10 2.0 3.8E-10 2.0 l.lE-09 2.0 



y— vel-gas pressure 

16 3.4E-09 - 7.6E-09 - 2.1E-08 - l.OE-08 - 2.3E-08 - 6.5E-08 

32 8.7E-10 2.0 1.9E-09 2.0 5.5E-09 2.0 2.7E-09 2.0 5.9E-09 2.0 1.7E-08 2.0 

64 2.2E-10 2.0 4.9E-10 2.0 1.4E-09 2.0 6.7E-10 2.0 1.5E-09 2.0 4.2E-09 2.0 

128 5.5E-11 2.0 1.2E-10 2.0 3.5E-10 2.0 1.7E-10 2.0 3.7E-10 2.0 l.lE-09 2.0 



y— magnetic 

16 7.0E-09 - 1.5E-08 - 4.3E-08 

32 1.8E-09 2.0 3.9E-09 2.0 l.lE-08 2.0 

64 4.4E-10 2.0 9.9E-10 2.0 2.8E-09 2.0 

128 l.lE-10 2.0 2.5E-10 2.0 7.0E-10 2.0 



Case B: 5 = 10-^ k = 27r(2, 1)/V5 
density-gas x— vel-gas 

16 7.7E-08 - 1.2E-07 - 2.4E-07 - 9.2E-08 - 1.4E-07 - 2.9E-07 - 

32 2.1E-08 1.8 3.4E-08 1.8 6.8E-08 1.8 2.6E-08 1.8 4.1E-08 1.8 8.2E-08 1.9 

64 5.6E-09 1.9 8.7E-09 1.9 1.7E-08 2.0 6.7E-09 1.9 l.lE-08 1.9 2.1E-08 1.9 

128 1.4E-09 2.0 2.2E-09 2.0 4.4E-09 2.0 1.7E-09 2.0 2.7E-09 2.0 5.4E-09 2.0 



y— vel-gas pressure 

16 2.3E-08 - 3.7E-08 - 7.8E-08 - l.OE-07 - 1.6E-07 - 3.5E-07 - 

32 5.6E-09 2.0 8.8E-09 2.1 1.7E-08 2.2 3.0E-08 1.8 4.7E-08 1.8 9.7E-08 1.9 

64 1.3E-09 2.1 2.1E-09 2.1 4.1E-09 2.1 7.8E-09 1.9 1.2E-08 1.9 2.5E-08 2.0 

128 3.2E-10 2.1 5.0E-10 2.1 l.OE-09 2.0 2.0E-09 2.0 3.1E-09 2.0 6.2E-09 2.0 



x— magnetic y— magnetic 

16 3.9E-08 - 6.3E-08 - 1.6E-07 - 4.7E-08 - 8.2E-08 - 1.9E-07 - 

32 9.0E-09 2.1 1.4E-08 2.2 3.3E-08 2.3 1.2E-08 2.0 1.9E-08 2.1 4.4E-08 2.1 

64 2.1E-09 2.1 3.3E-09 2.1 7.3E-09 2.2 2.9E-09 2.0 4.6E-09 2.1 l.OE-08 2.1 

128 5.2E-10 2.0 8.1E-10 2.0 1.7E-09 2.1 7.3E-10 2.0 l.lE-09 2.0 2.4E-09 2.1 
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Table 4: Convergence Rates Slow MHD Waves. 



Ncells Li Ri L2 R? ioc -^'o^ Li Ri L2 R2 -too Rc 



Case E:6 = 10-^ k = 27r(l,0) 
density-gas x— vel-gas 

16 7.5E-09 - 1.7E-08 - 4.7E-08 - 4.4E-09 - 9.9E-09 - 2.8E-08 

32 1.9E-09 2.0 4.2E-09 2.0 1.2E-08 2.0 l.lE-09 2.0 2.5E-09 2.0 7.0E-09 2.0 

64 4.8E-10 2.0 l.lE-09 2.0 3.0E-09 2.0 2.8E-10 2.0 6.2E-10 2.0 1.8E-09 2.0 

128 1.2E-10 2.0 2.7E-10 2.0 7.5E-10 2.0 7.0E-11 2.0 1.6E-10 2.0 4.4E-10 2.0 



y vel-gas pressure 

16 1.4E-08 - 3.1E-08 - 8.7E-08 - l.lE-08 - 2.3E-08 - 6.6E-08 

32 3.5E-09 2.0 7.7E-09 2.0 2.2E-08 2.0 2.7E-09 2.0 5.9E-09 2.0 1.7E-08 2.0 

64 8.7E-10 2.0 1.9E-09 2.0 5.5E-09 2.0 6.7E-10 2.0 1.5E-09 2.0 4.2E-09 2.0 

128 2.2E-10 2.0 4.8E-10 2.0 1.4E-09 2.0 1.7E-10 2.0 3.7E-10 2.0 l.lE-09 2.0 



y— magnetic 

16 l.lE-08 - 2.5E-08 - 7.1E-08 

32 2.9E-09 2.0 6.4E-09 2.0 1.8E-08 2.0 

64 7.2E-10 2.0 1.6E-09 2.0 4.5E-09 2.0 

128 1.8E-10 2.0 4.0E-10 2.0 l.lE-09 2.0 



Case F: 6 = 10-^ k = 27r(2, 1)/V5 
density— gas x— vel-gas 

16 6.1E-08 - 9.5E-08 - 1.9E-07 - 5.7E-08 - 8.9E-08 - 1.8E-07 - 

32 1.6E-08 2.0 2.4E-08 2.0 4.9E-08 2.0 1.4E-08 2.0 2.3E-08 2.0 4.7E-08 1.9 

64 3.9E-09 2.0 6.2E-09 2.0 1.2E-08 2.0 3.7E-09 2.0 5.8E-09 2.0 1.2E-08 2.0 

128 9.9E-10 2.0 1.6E-09 2.0 3.1E-09 2.0 9.2E-10 2.0 1.4E-09 2.0 2.9E-09 2.0 



y— vel-gas pressure 

16 1.4E-07 - 2.2E-07 - 4.4E-07 - 8.5E-08 - 1.3E-07 - 2.8E-07 - 

32 3.6E-08 2.0 5.7E-08 2.0 l.lE-07 1.9 2.2E-08 2.0 3.5E-08 2.0 7.0E-08 2.0 

64 9.1E-09 2.0 1.4E-08 2.0 2.9E-08 2.0 5.6E-09 2.0 8.7E-09 2.0 1.8E-08 2.0 

128 2.3E-09 2.0 3.6E-09 2.0 7.2E-09 2.0 1.4E-09 2.0 2.2E-09 2.0 4.4E-09 2.0 



x— magnetic y— magnetic 

16 6.1E-08 - 9.6E-08 - 2.0E-07 - 6.3E-08 - 9.9E-08 - 2.2E-07 - 

32 1.5E-08 2.0 2.4E-08 2.0 5.0E-08 2.0 1.6E-08 1.9 2.6E-08 1.9 5.5E-08 2.0 

64 3.9E-09 2.0 6.1E-09 2.0 1.2E-08 2.0 4.1E-09 2.0 6.5E-09 2.0 1.3E-08 2.0 

128 9.8E-10 2.0 1.5E-09 2.0 3.1E-09 2.0 l.OE-09 2.0 1.6E-09 2.0 3.3E-09 2.0 
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Fig. 9. — Loop advection: time evolution of the magnetic energy for the 2D (top) and 3D 
(bottom) cases, for three different resolutions. 




Fig. 10. — Loop advection in 3D: time evolution of B^, the magnetic field component aligned 
with the loop axis, for three different resolutions. 
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Table 5: Riemann Problem Set 

left-state right-state 



Test 


P 




Uy 




P 


By 




p 




Uy 




p 


By 


Bz 


Brio-Wu 


1 











1 


1 





0.128 











0.1 


-1 





Dai- Woodward 


1.08 


1.2 


0.01 


0.5 


0.95 


3.6 


2 


1 











1 


4 
\/47r 


2 


Ryu- Jones 


1 











1 








0.3 








1 


0.2 


1 





Fast Rarcf. 


1 


-2 








0.45 


0.5 





1 


2 








0.45 


0.5 






3.4. Orszag-Tang Vortex 



We now turn to a series of classical multidimensio nal tests for MHD co des. First, we 
compute the compressible Orszag-Tang vortex problem (lOrszag fc Tanglll979l ). We solve this 
problem on a computational domain of size L = 1 with periodic boundary conditions and a 
rectangular grid of 200x200 cells. We use an adiabatic index 7 = 5/3. We present results 
obtained using the PPM reconstruction scheme and primitive limiting, but very similar 
results are produced using characteristic limiting. The initial conditions in terms of the 
primitive variables are as follows 



W = [po, —uq sin(27r?/), Uo sin(27rx), 0, Pq, —Bq sin(27ry), Bq sin(47rx), 0]' 



(58) 



where po = 25/367r, Pq = 5/127r and Uq and B^ are defined in terms of the sonic Mach 
number, Ai = 1, and plasma beta, /3 = 10/3, respectively. Although the initial conditions 
are smooth, eventually the fiow develops a complex structure with sharp features and discon- 
tinuities. In Fig. [11] we show the contours of the numerical solution for density, gas pressure, 
kinetic energy and magnetic pressure at time t = 0.5. One dimensional cuts along the line 
y = 0.4277 for the thermal pressure (top) and magnetic pressure (bottom) are also shown in 
Fig. [12 The code maintains the symmetry of the solution with respect to the central point. 
In addition. Fig. [T2] shows that discontinuous features are captured within a few zones. Fi- 
nally, the solution ca n be compared with similar plots at the same solution time produced 



by other authors (e.g. iRyu et al.lll998t lT6thll2000l : IStone et al.ll2008l ). from which it appears 



that the code produces the correct fiow structures. 



3.5. Rotor Problem 



Another two-dimensional problem c ommonly used as a test for multidimensional MHD 
codes is the rotor problem described in iBalsara fc Spicerl (Il999l ). It consists of a rotating 
disk of dense material, threaded by a magnetic field initially directed along the x— axis, and 
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Fig. 11. — Orszag-Tang Vortex: thirty equally-spaced contour levels between max and min 
value of the numerical solution at t — 0.5 respectively for density (top-left), gas pressure 
(top-right), specific kinetic energy (bottom- left) and magnetic pressure (bottom- right). 




Fig. 12. — Orszag-Tang Vortex: One-dimensional cuts along the a;— axis at y = 0.4277 for 
the thermal pressure (top) and magnetic pressure (bottom). 
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embedded in a tenuous ambient medium at rest. As the rotor spins, it winds up the magnetic 
field fines, generating Alfven waves wfiicfi propagate into tfie ambient medium. Tfie problem 
is solved on a computational domain of size L = 1 witfi periodic boundary conditions and 
covered witfi a rectangular grid of 400x400 cells. Tfie adiabatic index is 7 = 1.4. We use 
tfie PPM reconstruction scfieme and primitive limiting, altfiougfi tfie use of cfiaracteristic 
limiting produces very similar r esults . Tfie s etup of tfie initial conditions corresponds to tfie 
first rotor problem discussed in iTotfil (l2000l ). i.e. 



W[x,t = 0] 



Wdisk if r < Tdisk, 

Wamb if r > Tamb, 



(59) 



wfiere W^i 



disk 



[pd^sk, -^{y-0.5), ^(x-0.5), 0, Po, Bo, 0, 0]^, Wamb = (po, 0, 0, 0, Pq, ^o, 0, 0)^, 
r = A/a;2 + y^, r^isk defines tfie disk radius and Vamb delimits tfie ambient medium. In 
tfie transition region between tfie ambient medium and tfie rotor tfie density and velocity 
fields are interpolated according to p = {pdisk — Po)f{f^) + Po, = —f{r)uo{y — 0.5) /r and 
= /(^)wo(a; - 0.5)/r, witfi /(r) = {vamb - r)/{ramb - Tdisk), wfiereas density and pressure 



remain uniform. Tfie results for tfiis test are presented in Fig. [131 wfiicfi effectively repro- 
duces Fig. 18 of iTotfil (120001 ). Tfie numerical solution appears very well befiaved and tfie 
code seem to pass tfiis test as well. 



4. Extension to Adaptive Mesh Refinement 



Following iBerger fc Colellal (Il989[ ) and iBalsaral (I2OOII ). we employ block-structured lo- 
cal refinement to increase the computational resolution where the accuracy of the solution 
ne eds to be improv e d. Ou r implementation is basically an extension of the MHD case 



of iMiniati fc ColeUal fl2007a[ ). 



Generalizing tfie case discussed in Sec. I2.1.H tfie problem domain is discretized on a 
fiierarcfiy of grids, r°...r^™"=, eacfi witfi its own spacing and refinement ratio nl^j- = 
h^/h^^^. We assume tfiat refinement ratio is always even. 



Calculations are performed on a fiierarcfiy of mesfies {f^ } 



f=0 



sucfi tfiat for eacfi 



Vt C r . Tfie base-level uniform rectangular mesfi spans tfie domain, so Vt^ 



CeHs 



for wfiicfi improved resolution is desired are marked for refinement, grouped togetfier into 
logically rectangular regions, and refined by a factor of n^^j to create tfie level 1 domain 
. Furtfier refinement levels, VL^, may tfien be created as needed in tfie same way starting 
from a refinement level VL^~^, witfi refinement ratio Tfie set of generated mesfies Vt^ 

is assumed to be properly nested, meaning tfiat (a) any control volume i G fi^ is eitfier 
completely covered by {nl^j)^ finer control volumes or by none, and (b) for any given pair 
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of meshes and there is always a layer of cells of separating the two. By analogy 
with the single grid case we can construct the set fi^^ corresponding to the faces in and 
likewise the set Q^'^'' for the corresponding edges. Similarly, for each level, we can define a 
divergence and curl operators for face and edge centered vectors, respectively. 

The part of an AMR level which is "covered" by refinement is denoted as the covered 
region, while the valid region of a given level i is the part of not covered by refinement. 
For computational convenience, solution values are maintained in covered regions as well as 
valid regions, but only the solution in valid regions is considered to be valid. The composite 
solution spans the computational domain and is the union of the valid-region solutions on 
each level. The coarse-fine interface between levels i and £ — 1 is denoted by dQ^. 



As in lBerger fc Colellal ( 119891 ). we refine in time as well as space, with At^~^^ = At^/n^^j. 
The update of the solution on the hierarchy of AMR levels can be described recursively as 
the update of a single AMR level i from time t^ to time t^ + At^ (Fig. H]). 

Extending the CT scheme described in this paper requires some additions to the stan- 
dard set of algorithmic tools ge nerally used for fu l ly cell -centered discretizations of hyperbolic 
conservation laws like that in iBerger fc Colellal ( 1l989l ). Most of the additional algorithmic 
pieces result from the addition of the solenoidal face-centered B field. 



4.1. Filling Ghost cells 

Before each single level update from time t^ to time t^ + At^, a ring of "ghost" cells 
sufficiently large to complete the stencils required to update valid-region data on each grid- 
patch is filled in. For the scheme described here, we require 6 and 4 ghost cells for PPM 
and PLM reconstruction methods, respectively. Where possible, ghost values are filled by 
copying valid-region data from other grids on the same level i or possibly by a discrete 
representation of physical domain boundary conditions. Where neither of these is possible, 
values must be interpolated from the coarser level £ — 1 solution. Cell-centered quantities are 
interpolated using a limited piecewise-linear scheme. The face-centered B field is likewise 
interpolated using a piecewise-linear scheme as follows: 

1. First, B^~^ is linearly interpolated in time to t^. (Recall that level i — 1 has already 
been updated from t^-^ to {t^-^ + At^-^), and t^'^ < t^ < {t^-^ + At^'^). 

2. Then, fine-level faces which overlie coarse-mesh faces are filled using piecewise linear 
interpolation of co-planar coarse-mesh faces. 

3. Finally, values for fine- level faces which do not overlie coarse- mesh faces are linearly 
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Fig. 13. — Rotor Problem: thirty equally-spaced contour levels between max and min value 
of the numerical solution at t = 0.15 respectively for density (top-left), gas pressure (top- 
right), Mach number (bottom-left) and magnetic pressure (bottom-right). 

LevelAdvance(£, At^ 

SingleLevelUpdate(£, t^, At'): 

Interpolate solution values as needed from next coarser level (i — 1) 

Update solution on level i using scheme described in Section 12.21 

if {i < ^max): increment fine flux registers 

if {i > 0): increment coarser-level flux registers 

t' := t' + At' 
Recursively update any finer levels: 

if {i < 4iax) then 
At^+i = 4*i 

re/ 

for n = 1, n^g^ 

Level Advance (^ + 1, At'+^) 
end for 

"Synchronize" levels £ and level £ + 1 
end if 
end LevelAdvance 



Fig. 14. — Recursive AMR timestep 
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interpolated between the two surrounding co-directional faces for which there are al- 
ready values (either fine-level faces on the coarse-fine boundary or fine-level values 
interpolated in step (1) ). 



Note that we hav e not found it necessary to use a divergence-free interpolation scheme like 
that described in iBalsaral ( 120011 ) to fill in values for ghost faces. 



4.2. Synchronization 



After the sub-cycled advance of level i + 1, the solutions on levels i and i + 1 have 
reached the same solution time (t^ = t^^^), and are then syn chronized. For ce l l-cent ered 
conserved variables, synchronization is identical to that used in lBerger fc Colellal (Il989[ ): 



1. Replace level i solution with the averaged level i + 1 solution in covered regions. 

2. Because the fluxes used to update the fine-level solution were computed independently 
of those used to compute coarse-level updates, con serva tion will not be rnaintain ed at 
coarse- fine interfaces. In iBerger fc Colellal ( 119891 ) and Martin fc Colellal (l2000l ). flux 
registers are defined along the coarse-fine interface between levels i and £ + 1, in which 
the fluxes used to compute coarse- and fine-level updates are stored. Since we consider 
the fine-level fiuxes to be more accurate, we update the solution in the coarse cells 
adjoining the coarse-fine interface with the "refiux-divergence" of the difference of the 
fluxes: 



(60) 



n=l 



where is the vector of conserved quantities, Dji is the "reflux divergence" operator, 
is the vector of fluxes used to update U^, is the average of the level (£ + 1) 

fluxes on the underlying level i faces, and the sum is over sub-cycled i + 1 timesteps. 



Synchronization for the f ace-ce ntered magnetic fleld looks similar, and takes the same 
form as described in IBalsaral (120011 ) : 



1. Replace level i magnetic fleld with the averaged level £ + 1 solution on covered faces. 

2. Because the edge-centered electric flelds on each level are computed independently, 
the composite B fleld will no longer be solenoidal. We treat this using an analogue 
to the face-centered flux registers used for cell-centered data, as presented in iBalsara 
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(l200ll ). We store the edge-centered electric fields along coarse- fine interfaces, then 
increment the coarse-level magnetic field at faces bordering the coarse-fine interface 
with a reflux-curl operator applied to the coarse-fine mismatch in the edge-centered 
electric fields. 



4.3. Regridding 

It is often desirable for refined regions to periodically adapt as the solution evolves in 
time. When newly refined regions are created, cell-centered fields are also interpolated from 
coarse-level data using limited piecewise- linear interpolation. For B field values o f newl y 



refined faces, we use the divergence-free interpolation scheme described in iBalsaral ( 1200 ll ). 
The scheme is defined for refinement ratios Uref = 2, but it can be applied recursively for 
larger values of the refinement ratio. 



4.4. Shock-Cloud Interaction 



As an example of an AMR-MHD application we compute the interaction of a cloud with 
a strong shock wave. This process is common in the interstellar medium where s hocks pro- 



duced by supernova explosio ns interact with the surrounding multi-phase medium ( iKlein et al. 



I994J : iMac Low et al.l Il994l ). Related processes, characterized by similar hydrodynamic 



medium (Schiano et al. 1995 : Jones et al. 19961: Vietri et al.ll997j : 



1997 



igggl - booohT or supersonic clouds collisions (iLattanzio et al 



1985 



1999aD. 



Miniati et al.lll999bl: 


Greeori et al. 


Klein et al. 


I995I; 


Miniati et al.l 



The initial conditions for our problem are as follows: the background gas has unit density 
and thermal pressure, and is at rest; a cloud with the same pressure but 10 times higher 
density, moves through the thin gas with a velocity Vc = —3.47871373 along the x direction. 
A plane-parallel shock with Mach number = 10 propagates along the same axis but in the 
opposite direction to the cloud. The initial magnetic field is uniform, of unit strength and 
aligned with the x— axis. The computational box has dimensions [0, 1] x [0, 0.5] x [0, 0.5]. The 
boundary conditions correspond to supersonic inflow and outflow for the lower and upper 
boundaries of the x— axis, and are periodic otherwise. The calculation is carried out in three 
dimensions on a base grid of 64 x 32 x 32 cells. Two additional levels of reflnement, with 
reflnement ratio 4, are generated dynamically in regions where the normalized undivided 
density gradient \Ap\/p > 0.1, and/or in the presence of shocks according to the criteria 
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|AP|/P > 0.1 and V ■ V < 0. 

Fig. [TSl shows from top to bottom the solution for the density, gas pressure, magnetic field 
magnitude and plasma beta parameter (/3 = Pgas/Ps), at time t = 0.021251, corresponding 
to 160 cycles on the finest level. The main features, discussed at length in the above papers, 
are correctly reproduced. The plane shock front moving from the left crushes the cloud. As 
the cloud moves to the left, it creates a bow shock in front of it, where the pressure has its 
highest value. The cloud motion also generates a low pressure region at its rear, where the 
magnetic pressure dominates the gas pressure and the beta plasma is much less then unity. 
In addition, as the fluid flows past the cloud, the magnetic field line s entrained in the cloud 



body fold on themselves creating a current sheet (I Jones et al.lll996l ) 



Note that the maximum value of the normalized divergence of the magnetic field | AxV ■ 
B|/|i?|, is a fewxlO"^^. This is completely negligible with respect to the solution value and 
demonstrates that our implementation of the above operators for the coarse fine magnetic 
field interpolation and refiuxing operations is correct. 



5. Extension to Cosmology 

We now describe the extension of our MHD algorithm to the case of cosmological appli- 
cations. This will in clude only a basic descrip tion of the cosmological code, for it is presented 



in detail elsewhere (IMiniati fc Colellall2007al ). 



For cosmological simulations it is preferable to transform away the expansion of the uni- 
verse through the use of a comoving frame of reference. Thus we operate the transformation 

X ^ a(ty^ X (61) 

from lab to comoving coordinates, where a{t) is a function of time that defines the physical 
size of spatial scales, and a/ a is the expansion rate of the universe. In addition we subtract 
out the velocity component due to the expansion of the universe, and retain the peculiar 
proper velocity, i.e. 

u ^ u — ax. (62) 

Finally, it is convenient to use comoving density and pressure, i.e. those expressed in terms 
of the comoving volume x^ as opposed to the proper volume a^a;^. 



p ^ a^ p, 
p a^p. 



(63) 
(64) 
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Fig. 15. — Shock-cloud interaction: from top to bottom, solution for density, gas pressure, 
magnetic field magnitude and plasma beta parameter, at time t = 0.021251. 
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Similarly, although the magnetic flux naturally scales with a^(t), we transform it to a pseudo- 
comoving variable 



B B. 



(65) 



The above transformations allow writing the conservation and induction equations in a form 
that, except for the appearance of source terms, closely resembles the original ones. This 
similarity allows us not only to solve for the MHD system of equations in the cosmological 
framework with the same algorithm described in Sec. [2] but also to apply the same menagerie 
of AMR tools described in Sec. IHwith virtually no modiflcation. 



In fact, in the comoving frame, x, the conservation equations read 



dU 1 _ ^ 

at a{t) 



(66) 



where U and F are deflned exactly as in (fTT!) but are now expressed in terms of peculiar 
velocity, comoving density, comoving pressure and pseudo-comoving magnetic field. The 
source term on the RHS is 



S{U) 



a 
a 





puo 



\ 



2pe-f 



+ p 



f \ 

9D-1 
u • g 



V / 



(67) 



where the first term, oc a/a, accounts for adiabatic losses of momentum, energy, and magnetic 
field, and the second term, oc g, is due to gravity. Similarly, we rewrite Faraday's law in the 
comoving frame in terms of the peculiar velocity and pseudo-comoving magnetic field, i.e. 



dt 



Id la 

{UjBd - BjUd) - 2~Bd. 



a dxn 



(68) 



Based on Eq. fl66l) and (|68|) the time update of U and B is then done according to 

1 At 



U. 



n+l 



B 



n+l 



1 

2 



(V ■ F)T-^ 



in+l 



1 At 



('a''+3a"+i)5 Ax 



[D X Ef^'^ 



d,i+\e<i ' 



(69) 
(70) 
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where the time centered fluxes and electric fields, as well as the synchronization between 
face and cell centered magnetic field variables, are computed using the algorithm defined in 
Sec. [2J A second order estimate of the source term can be obtained by the simple time average 
5"+2 ~ 1(5" + ^^^+1). In reality, the source terms associated with gravity and expansion are 
implemented using a slightly more sophisticated method that estimates the change in kinetic 
energy due to the work by gr avity, directly frorn . the cha nge in the momentum components. 
This is described in detail in iMiniati fc Colellal (l2007aj). Similarly, after the face-centered 
magnetic field variables have been updated in time, the cell-centered values are synchronized, 
and the change in magnetic field energy due to cosmic expansion is computed from the 
corresponding change in the magnetic field components. 



5.1. MHD Santa Barbara Test 

In this final test we present an MHD version of the 'Santa Barbara Cluster Comparison 
Project'. The tests consists of the formation of a massive cluster of galaxies in a 64 Mpc 
volume. The cosmological model is an Einstein-De Sitter universe (i.e. with critical total 
matter density) with 10% of the total matter density in baryons, and an ex pansion rate 
given by a Hubble parameter, Hq = 50 km s~^ Mpc~^ (see additional details in iFrenk et al. 



( 119991 )). The purpose of this calculation is to test our MHD solver in a realistic cosmological 
application. To compare with previously published results, the dynamic role of the magnetic 
field remains negligible throughout the calculation, which we ensure by adopting a sufficiently 
small initial magnetic seed. The geometry of such fields is immaterial, and is chosen to be 
uniform for convenience. 



The calculation is performed basically as in iMiniati fc Colellal (j2007al ). except for the 
initial redshift which is z = 30 (instead of 40). So, two initial grids are in place at simulation 
start: a base grid covering the entire 64 Mpc^ domain with 64^ cells and 64^ particles; and 
a second grid, also with 64^ cells and 64^ particles, but only 32 Mpc on a side and placed in 
the central region of the base grid, thus yielding an initial cell size of 0.5 Mpc. Refinement is 
applied only within the latter higher resolution region to cells with a total mass larger than 
6.4 X 10^° Mq. We allowed for 5 additional levels of refinement (for a total 6 level hierarchy), 
with a constant refinement ratio riref = 2. The size of the finest mesh is about 15 comoving 
kpc. The timestep is limited by the most stringent among the following three conditions: 
the Courant-Friedrichs-Lewy condition on the MHD waves, with coefficient Ccfl = 0.8, an 
analogous CFL condition based on the speed of the collision-less particles, with coefficient 
Cpart = 0.5, and the requirement that the expansion of the universe during a time-step is 
less than 2%. The calculation was performed using the PPM reconstruction scheme and the 
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HLLD Riemann solver. 



The resu lts of the calculat i on are summarized by the radial plots in Fig. [ini where 
in analogy to iMiniati &: Colellal (l2007aj) we show results from two other simulation codes. 
As expected, the MHD solver (in the lir nit of vanishing magnetic field) produces virtually 
the same results as the hydro solver in iMiniati fc Colellal (j2007al ). attesting therefore to 
the same reliability for highly nonlinear calculations, involving high Mach number flows 
and large dynamic range in the fluid quantities. Finally, the left panel of Fig. [T7] shows 
the magnetic field magnitude (in arbitrary units) distribution on a slice passing through 
the cluster center. The magnetic field is stronger in the core region where it also shows 
substantial spacial structure down to the smallest resolvable scales. On the left panel of 
the same Figure we present the distribution of the normalized magnetic field divergence, 
|(Ax/-B)V-B|. The bulk of the distribution is at the level of 10~^^ or so, and the max value 
is around 10^^^. Again, this value is completely negligible with respect to the solution value. 
While larger than the value obtained in the previous test example, this is expected given the 
much larger number of integration steps (about 10^) in the current case. 



6. Summary 



We have presented the implementation of a three-dimensional scheme for MHD in the 
AMR code CHARM. The scheme uses a hybrid discretization, in the sense that fluid quantities 
are cell-centered, and magnetic field variables are face-centered. The algorit hm is based o n 
the full 12-solve spatially unsplit Corner- Transport-Upwind (CTU) scheme (jColellalll990l ). 
The fluid quantities are updated using the PPM method, while the magnetic field evolution 
is computed through a CT r nethod. The edge - center ed electric fields necessary for the 
CT step are computed as in iGardiner fc Stond (120051 ). We employ a simplified version 
of the m ultidimensional MHD source t erms required i n the predictor step for high-order 
accuracy (IGardiner fc Stond l2005l . 120081 ) . which is as in lCrockett et al.l (|2005[ ). This greatly 
simplifies the three-dimensional version of the algorithm with respect to the original form, 
without compromising the accuracy and robustness of the solutions. 

The algorithm is implemented in an AMR framework. This requires synchronization 
operations across refinement levels, including face-centered restriction and prolongation oper- 



ations and a reflux-cur! 



oper ation, which is necessary to maintain a divergence-free magnetic 



field solution iBalsaral (1200 ll ). The code works with any even refinement ratio, although val- 



ues above 4 are unusual. Our tests demonstrate that the code converges at the expected 
rate, is robust in problems involving strong shocks, maintains the magnetic field divergence 
at a negligible value and is suitable for astrophysical and cosmological applications. 
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Fig. 16. — Radial profile of dark matter (top left), baryonic gas (top right) temperature 
(middle left), baryonic fraction (middle right), radial velocity for dark matter (bottom left) 
and gas (bottom right). In addition to the results from CHARM ( open circles), for compari son 
we also show those from the ENZO AMR code f filled triangles') [Bryan fc Norman! (119971 ) as 
well as those from the HYDRA SPH code (open stars) Couchman et al. f 19951 ). 
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Fig. 17. — Left: Distribution of the magnetic field magnitude on a plane across the center of 
the simulated cluster. Right: Histogram of the absolute value of the normalized divergence 
of the magnetic field, |(Ax/-B)V-B|. 



